Monte Carlo calculations for liquid 4 He at negative pressure 
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Abstract 



A Quadratic Diffusion Monte Carlo method has been used to obtain the 
equation of state of liquid 4 He including the negative pressure region down to 
the spinodal point. The atomic interaction used is a renewed version (HFD- 
B(HE)) of the Aziz potential, which reproduces quite accurately the features 
of the experimental equation of state. The spinodal pressure has been calcu- 
lated and the behavior of the sound velociy around the spinodal density has 

been analyzed. 
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Recent experiments [|T]-0] on cavitation in liquid 4 He at low temperatures have motivated 
the theoretical study ||-[J of liquid helium properties at low temperature and negative pres- 
sure. Some interesting questions have thus been raised, as the determination of the tensile 
strength (i.e., the magnitude of the negative pressure at which cavitation becomes likely), 
and the spinodal pressure (i.e., the pressure at which liquid helium becomes macroscopi- 
cally unstable against density fluctuations). The standard (or classical) theory of nucleation 
predicts for liquid helium [|l(| a tensile strength rising from 8 atm at 2 K to 15 atm at 0.5 
K, and the experiment around 1.5 K reported in Ref. seemed to confirm this prediction. 
However, Maris and Xiong estimated 0|| the spinodal pressure to be -9 atm at K, so 
that the predictions of the standard theory must be incorrect (since the liquid cannot exists 
for pressures lower than the spinodal one). They also carried out an experiment 0, whose 
results are in contradiction with those of Ref. Jl|, obtaining lower values for the tensile 
strength. More recently, several experiments providing information about the cavitation 
problem have been performed. Unfortunately, as they do not rely on an accurate pressure 
calibration, no tensile strength values are reported. 



In Refs. jm the spinodal pressure was estimated by fitting to the measured ]TJJ sound 
velocities c as a function of pressure P several polynomial and Pade forms, and then extrapo- 
lating into the negative pressure region to determine the zero of c(P). From a different point 
of view, the spinodal pressure was calculated in Ref. using two different phenomenolog- 
ical models that reproduce the equation of state in the measured positive pressure region. 
Although an overall agreement between the phenomeno logical calculations and the empirical 
results was obtained, some questions arise, as for instance to what extent the extrapolated 
results depend on the form used in the fit, or on the density functional used in the calcula- 
tions. It is therefore necessary to handle with a precise equation of state for liquid 4 He valid 
in the full range of pressure, down to the spinodal one. 

Many-body techniques have achieved a high level of accuracy in the description of liquid 
helium. In particular, Monte Carlo (MC) methods JT2] give exact information, apart from 



statistical uncertainties, on the ground state of bosonic systems both at zero and finite 



temperature. The well known interatomic interaction for helium has been a key ingredient 
to reach an excellent agreement between the theoretical results and the experimental data. 
Recently, we have used a Quadratic Diffusion Monte Carlo (QDMC) method to calculate 
the equation of state in the positive pressure region . One of the main conclusions of Ref . 
TB| is that the HFD-B(HE) potential suggested by Aziz et al. |14}| , hereafter referred to as 



Aziz II potential, improves the results obtained with the Aziz potential ||15|| , especially when 
the density dependence of derivative magnitudes of the energy is considered. In this work, 
we have extended those calculations to lower densities, with the hope of studying without 
ambiguity the zero-temperature properties of homogeneous liquid 4 He in this zone. 

The QDMC method solves stochastically the Schrodinger equation in imaginary time 



assuming a short-time approximation form for the Green's function [|T6|. The ground-state 
wave function is sampled in an iterative process after a time larger enough to project out 
higher energy components. Rigorously, the exact ground-state energy is obtained when the 
limit At — > is considered. In linear DMC algorithms one has to perform calculations at 
several time steps, and then extrapolate to the exact value. The QDMC method, which 
has evidenced a quadratic dependence on At [TB|,[IjJ , improves the efficiency of the diffusion 
algorithm making feasible to use larger time steps than in DMC and avoiding the necessity 
of the extrapolation to At = 0. 



In order to guide the diffusion process a Jastrow trial wave function of the form |T8 

* t (r) = n ex p 
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has been used in our calculations. The values of the parameters appearing in Eq. ([[]) 
(L = 0.2, A = 2.0 a, A = 0.6 cr, b = 1.20 a , where a = 2.556 A) have been fixed to optimize 
the Variational Monte Carlo (VMC) estimation of the energy at the equilibrium density 
(po = 0.365 cr -3 ). The Monte Carlo simulation has been carried out with 108 particles for 
p < 0.328 cr -3 and 128 particles for p > 0.328 a" 3 , maintaining an asymptotic population 
of 400 walkers. The errors associated to the use of both a finite volume simulation box and 
a finite population have been analyzed and, in all smaller than the size of the 



inherent statistical fluctuations. 

In Table I are reported the total energies (in K) obtained with the Aziz II potential, for 
different values of the density (in units of cr~ 3 ). The experimental values of Ref. [T{| are also 



displayed. The origin of the slight differences observed between theory and experiment was 
discussed in Ref. [T3[ . 



Derived quantities of the energy such as the pressure or the sound velocity have been 
obtained through a third and fourth degree interpolation, with unnoticeable changes when 
larger degrees were introduced in the interpolation method. The QDMC prediction of P(p) 
is shown in Fig. 1 (solid line) for the whole range of densities, in comparison with experi- 



mental data for positive pressures [0. The agreement between the Aziz II results and the 
experimental data is quite impressive. 

At this point, we would like to draw the attention on the quality of the extrapola- 
tions coming from the previously available data, laying mainly in the positive pressure 
region. In the majority of microscopic calculations on liquid helium the energy per particle 
is parametrized using a polynomial form 

E/N = e + B(^)\ C ( P -^) 3 . (2) 

On the other hand, in calculations based on Density Functional Theory, the form 

E/N = bp + cp 1+1 (3) 



proposed by Stringari p0| , has proved to be very efficient in describing properties of ho- 
mogeneous and inhomogeneous (including an additional surface term) liquid 4 He. We have 
used both forms to fit our previous QDMC results ||13|| , which included only one point below 
the equilibrium density. Proceeding in this way (i.e., considering only the five last densities 
of Table I), we have found that although both fits are compatible with the previous results 
of the energy, only the second form, see Eq. (0), predicts the new QDMC results at densities 
lower than 0.328 a~ 3 . This fact is reflected in Fig. 1, where the pressure derived from these 
fits is plotted as a function of density. The short-dashed line corresponds to the fit (|2|) and 



the long-dashed line to the fit (|3|). The starting points from the left of all curves depicted in 
Fig. 1 correspond to the location of the spinodal point. The QDMC result for the spinodal 
pressure is P c = — 9.30±0.15 atm, corresponding to a density p c = 0.264±0.002 a~ 3 . The fit 
(|) predicts (p c , P c ) = (0.266, -9.08) whereas the fit (g) gives (p c , P c ) = (0.292, -6.66). This 
Figure illustrates that the extrapolation from the positive pressure region to the negative one 
is quite sensitive to the form used in the fit. One can also see the capability of the built-in 
density functional given in Eq. (|]) to reproduce quite accurately the equation of state of 
liquid 4 He at any density, even when only positive pressure values are used as input in the 
numerical fit. It is also noticeable the stability of Eq. (§) when all the QDMC energy results 
(Table I) are used to fix the optimal parameters b, c and 7, being their relative changes less 
than 5 %. In this case, the spinodal point turns out to be (p c , P c ) = (0.267, —9.16). If the fit 
(0) is extended to all the energy values, the spinodal point (p c , P c ) = (0.267, —9.24) comes 
close to the Monte Carlo prediction. 

In Figure 2 is displayed the sound velocity c as a function of pressure P. The points are 



the experimental values of Ref. [II], and the solid line corresponds to the QDMC results. 
The accuracy provided by the Aziz II potential is again remarkable, giving results for the 
sound velocity in close agreement with the experiment. It can be seen that c drops to zero 
very fast when approaching the spinodal point. The behavior of c near P c is expected to 
be of the form c oc (P — P C Y, being v the critical exponent. It is known |7] that v = 1/4 
provided that the quantity 

Hf). 

be different from zero. Our QDMC estimation of P c " is 1200 ± 100 Ka 3 . Therefore, the clear 
departure from zero of P c " guarantees that v — 1/4, in disagreement with the Maris' model 
H which, explicitly taking the value P c " = 0, predicts v = 1/3. Nevertheless, in the positive 
pressure region the experimental values are very well reproduced by the above form with v 
close to 1/3. In fact, we have found that for pressures only slightly higher than P c and up 
to almost solidification pressure, the behavior is also c oc [P — P c ) u with the exponent given 



by 

V = i/ = -P c K Mo , (5) 



where is the isothermal compressibility 

] 

p \dP 



and Uq is the Griineisen constant 



P dc 

u=- — , 7 
c ap 



both quantities evaluated at the saturation density p . Our QDMC results give for these 
quantities the values k = 0.01231 ± 0.00005 atm" 1 and u = 2.81 ±0.04 with p = 0.3645 ± 



0.0003 a 3 , to be compared with the experimental values [11] k^ p = 0.01230 atm 1 , Mq XP 



2.84 and po* P = 0.3646 a" 3 . The QDMC result for the exponent (|) is u = 0.322 ± 0.007, 
close to the value 1/3 pointed out in Refs. 

The power-law behavior of c is shown in Figures 3 and 4. The quantity (c/cq) 1 ^, where 
Co is the sound velocity at the equilibirium density, is displayed in Figure 3 against 1 —P/P c . 
When v = uq a straigth line is obtained in a very wide range of pressures (Fig. 3), whereas 
the case v — 1/4 manifests a nearly quadratic behavior. However, for pressures close to the 
spinodal pressure one can see (Fig. 4) that the linear behavior is obtained taking the critical 
exponent as 1/4. 

To summarize, we have calculated the equation of state of liquid 4 He in a QDMC frame- 
work using the Aziz II potential from the spinodal point up to the solidification point. The 
characteristics of the spinodal point have been evaluated without the uncertainty of the ex- 
trapolation from the positive pressure data. As a byproduct, we have noted that our QDMC 
results are very well fitted by the form (|3|), giving thus additional support to this purely 
phenomenological density functional. 

This work has been partly supported by DGICyT (Spain) Grant Nos. PB90-06131, 
PB92-0761 and PB92-082. Most of the simulations were performed on a multiprocessor 
CM2 of the CEPBA (Centre Europeu de Paral-lelisme de Barcelona). 
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FIGURES 

FIG. 1. Pressure as a function of density. Points: experimental results solid line: QDMC 
result, short-dashed line: fit (||), and long-dashed line: fit (|3|). The fits (y) and @ include only 
the results from .328 to .438 cr~ 3 . 

FIG. 2. The sound velocity as a function of pressure. The experimental points are taken from 



Ref. [11 1, and the solid line corresponds to the QDMC results. 



FIG. 3. Power-law behavior of the sound velocity using as a critical exponent vq @ and 1/4. 

FIG. 4. Same as Figure 3 for pressures close to the spinodal one. To illustrate more clearly the 
power-law behavior, the quantity (c/cq) 1 ^ / (1 — P/P c ) has been plotted against 1 — P/P c for the 
values v = vq and 1/4. 



9 



TABLES 

TABLE I. Results for the total energies (in K) at several densities (in <r~ 3 ) from the QDMC 
calculations with the Aziz II potential and experiment IE]. 



p E/N E/N (exp) 

0.264 -6.376 ± 0.010 

0.267 -6.441 ± 0.010 

0.285 -6.692 ± 0.010 

0.300 -6.892 ± 0.010 

0.310 -7.011 ± 0.010 

0.328 -7.150 ± 0.010 

0.365 -7.267 ± 0.013 -7.17 

0.401 -7.150 ± 0.016 -7.03 

0.424 -6.877 ± 0.022 -6.77 

0.438 -6.660 ± 0.017 -6.55 
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